/* program moulton3.sas -- this program computes intrclass correlations

   jda, 5-5-99
   THE FOLLOWING MACRO VARIABLES MUST BE PASSED
   FROM THE "HOST" PROGRAM;

%let randvar = ; * defines groups within which errors are correlated;
%let k = ; * number of regressors (including intercep);
%let weight = ; * dummy for whether regressions are weighted;     
%let flag = ; * to signal residual calculation;                   */

data residvar;
  set test2;
  sum_sq = v_ij*v_ij;
  keep sum_sq;
  run;

proc summary data=residvar;
  var sum_sq;
  output out=residvar sum=;
  run;

data residvar; 
  set residvar;
  var_v = sum_sq/(_freq_-&k);
  keep var_v;
  run;

proc summary data=test2 nway;
  class &randvar;
  var v_ij;
  output out=sigma1 sum=sum_j;
  run;

data test2;
  merge test2 sigma1;
  by &randvar;
  sum_k = v_ij*(sum_j-v_ij);
  drop _type_ _freq_ sum_j;
  run;

proc summary data=test2 nway;
  class &randvar;
  var sum_k;
  output out=sigma2 sum=sum_jk;
  run;

* sum_jk is the double sum over j and k (j ne k) of (v_ij)*(v_ik);
* ignore v_bar in Moulton formula since it's zero by construction;

data sigma2;
  set sigma2;
  denom = _freq_*(_freq_ - 1);
  dummy = 1;
  drop _type_;
  run;

* now sum over groups (i) to get numerator and denominator of rho_hat;

proc summary data=sigma2 nway;
  var sum_jk denom;
  output out=sigma3 sum=num denom;
  run;

data ds&step;
  merge sigma3 residvar;
  rho&step = num/(var_v*denom);
  dummy = 1;
  keep rho&step var_v dummy;
  run;

proc print data=ds&step;
  var rho&step;
  run;






